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Two  new  algorithms,  the  Dual  Active  Set  Algorithm  for  linear  programming  and 
the  Conjugate  Gradient  Active  Set  method,  will  be  presented  for  solving  linear  and 
quadratic  programming  problems  respectively.  Both  methods  use  an  active  set  strat- 
egy. The  first  method  identifies  the  active  indices  by  maximizing  an  unconstrained 
dual  function  in  a finite  number  of  iterations,  and  the  second  method  solves  the  primal 
problem  in  a finite  number  of  iterations.  We  implement  the  dual  active  set  algorithm 
on  some  test  problems  extracted  from  the  Netlib  collection  of  linear  programming 
test  problems.  We  compare  the  relative  speed  of  some  gradient-based  algorithms  and 
the  conjugate  gradient  active  set  algorithm  using  a test  problem  from  engineering. 
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CHAPTER  1 
INTRODUCTION 

1.1  Overview 

The  notion  of  the  active  set  strategy  has  been  widely  embedded  in  a vari- 
ety of  algorithms  for  solving  various  optimization  problems.  It  is  frequently  used 
as  the  framework  in  many  algorithms  for  quadratic  programming  with  inequality 
constraints. 

For  a minimization  problem  with  inequality  constraints,  the  active  set  strategy 
divides  the  inequality  constraints  into  two  sets  and  then  solves  a sequence  of  sub- 
problems with  a set  of  equality  constraints  throughout  the  iterations.  During  the 
iterations,  the  procedures  for  adding  and  dropping  equality  constraints  to  these  sub- 
problems are  adopted  to  minimize  the  objective  function  of  the  problem.  Most  of  the 
algorithms  employing  the  active  set  strategy  are  designed  to  solve  the  primal  mini- 
mization problem.  A brief  literature  review  for  them  can  be  seen  in  section  (1.2).  To 
distinguish  the  active  set  strategy  used  for  solving  the  primal  minimization  problem 
from  another  type  of  active  set  strategy  which  solves  the  dual  problem,  we  call  it  the 
primal  active  set  strategy  in  the  remainder  of  this  paper. 

Recently  a new  active  set  scheme,  the  Dual  Active  Set  Algorithm,  has  been  de- 
veloped by  William  W.  Hager  [15].  Instead  of  solving  the  primal  problem  as  in  the 
Primal  Active  Set  Algorithm,  the  dual  active  set  algorithm  solves  the  Lagrangian 
dual  problem  of  the  given  problem.  Contrary  to  adding  in  constraints  for  the  Primal 
Active  Set  Algorithm,  the  dual  active  set  algorithm  drops  constraints  during  the  it- 
erations. This  new  algorithm  brings  in  another  direction  in  the  research  of  active  set 
strategy. 
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The  primary  goal  of  this  paper  is  to  introduce  two  new  algorithms,  the  dual  active 
set  algorithm  for  linear  programming  and  the  conjugate  gradient  active  set  algorithm 
for  quadratic  programming.  These  two  algorithms  employ  the  above  two  active  set 
strategies  respectively.  The  first  algorithm  modifies  the  dual  active  set  algorithm  by 
Hager  [15]  to  solve  a linear  programming  problem,  the  latter  one  adopts  the  primal 
active  set  strategy  to  minimize  the  quadratic  subproblem  along  the  conjugate  gradient 
direction  while  maintaining  a set  of  constraints.  Numerical  experiments  are  provided 
for  both  methods. 

For  the  remaining  two  sections  of  this  chapter,  we  describe  the  concepts  of  the 
Primal  Active  Set  Algorithm  and  the  Dual  Active  Set  Algorithm.  In  chapter  2, 
we  introduce  the  dual  active  set  algorithm  for  linear  programming.  Two  versions 
of  the  method  are  provided  and  discussed  followed  by  theorems  for  algorithms  to 
stop  in  finitely  many  iterations.  Numerical  experiments  on  problems  from  Netlib  are 
also  included  at  the  end  of  this  chapter.  We  present  the  conjugate  gradient  active 
set  algorithm  for  quadratic  programming  in  chapter  3.  Chapter  4 contains  a brief 
conclusion. 

1.2  Primal  Active  Set  Algorithm 

The  idea  of  employing  the  active  set  strategy  in  an  algorithm  to  solve  a min- 
imization problem  with  inequality  constraints  is  described  in  the  following.  Starting 
from  a feasible  point  of  the  feasible  domain  defined  by  the  inequality  constraints, 
the  algorithm  identifies  a set  of  active  constraints  called  the  working  index  set.  In 
the  successive  sub-iterations,  the  algorithm  minimizes  the  problem  on  the  subspace 
defined  by  the  constraints  with  indices  in  the  working  index  set.  During  the  line 
search,  new  constraints  may  be  added  into  the  current  working  index  set  whenever 
some  non-active  constraints  are  reached  by  the  search  direction.  When  the  algorithm 
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attains  a minimizer  of  the  subproblem  with  the  constraints  corresponding  to  the  cur- 
rent working  index  set,  certain  optimal  conditions  are  employed  to  drop  indices  from 
the  working  index  set  if  that  minimizer  does  not  satisfy  the  optimal  conditions  for 
the  original  problem. 

We  state  some  common  characteristics  of  those  algorithms  employing  the  primal 
active  set  strategy  for  solving  a minimization  problem  with  inequality  constraints  as 
follows: 

1.  Subproblems  with  different  binding  constraints  are  solved  through  the  iterations 
to  decrease  the  objective  function. 

2.  Each  point  generated  by  the  algorithm  is  feasible  in  the  subspace  defined  by 
the  constraints. 

3.  Finite  termination  of  the  algorithm  can  be  shown  if,  for  each  corresponding 
active  working  set,  the  strict  descent  of  the  objective  function  at  each  main 
iteration  is  assured. 


To  expose  the  idea  of  the  primal  active  set  strategy,  we  present  a primal  active  set 
algorithm  for  a specific  type  of  minimization  problem.  This  algorithm  is  an  extended 
version  of  the  primal  active  set  algorithm  given  in  Hager  and  Shih  [18]. 

Consider  the  following  problem 

min  f(x) 


s.t.  I < x < u 


(1.1) 


where  / is  a convex  function,  / and  u are  constant  vectors. 

In  the  following  algorithm,  the  sub-iteration  between  the  main  iterations  xjt  and 


£;t+1  is  denoted  by  yj. 
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Primal  Active  Set  Algorithm 

Initialization:  Let  x0  be  feasible  for  the  problem  (1.1),  set  j = 0,  set  y0  = Xk,  and 
define 

Ibq  = { i : (xj,),'  — or  ( Xk)i  = u,’  }. 

Sub-iteration  : Let  x satisfy  either  (a)  or  (b)  below 

(а)  If  x,  < /,  or  x,  > Ui  for  some  i in  the  complement  of  Is,,  then 
/(x)  < f(yj)  when  j = 0,  and  /(x)  < f(yj)  when  j > 0. 

(б)  If  / < x < u,  then  x = arg  min{/(x)}  subject  to  the  active  constraints 
corresponding  to  Ib}- 

If  (a)  holds,  then  set  = yj  + a(x  — yj)  where  a is  the  largest  a such  that 

/ < yj  + a(x  — yj)  < u for  all  0 < a < a,  increment  j,  and  set 


hi  = { * : ( Vj)i  = h or  (y^,-  = u{  }, 


proceed  to  the  next  sub-iteration. 

If  (6)  holds  and  Sjf(x)  satisfies  (1)  and  (2) 

(1)  > 0 for  i € Ib,  and  (x.)  = 

(2)  < o for  i € Ib,  and  (x.)  = ux. 


then  x is  an  optimal  solution  to  (1.1). 

Otherwise,  let  denote  the  final  set  Is,  generated  in  the  sub-iteration, 
set  I B0  = where  % G I^k  is  chosen  to  satisfy  the  following  conditions: 


dyf(x) 

dxi 

dy  fix) 

dxi 


> 0 if  X,-  = Ui 
< 0 if  x,  = U 
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set  Xfc+i  = x , increment  k,  and  return  to  the  sub-iteration. 


The  Simplex  method,  developed  by  George  B.  Dantzig  in  1947  for  the  linear  pro- 
gramming problem,  is  essentially  a special  type  of  active  set  method.  It  activates 
a constraint  and  deactivates  another  at  each  iteration  after  a vertex  of  the  feasible 
domain  is  found  by  a Phase  I procedure.  A non-simplex  active-set  method  for  linear 
programming  can  be  found  in  Gill,  Murray,  and  Wright  [10].  Many  algorithms  which 
employ  the  active  set  strategy  have  appeared  in  literature.  For  example,  Mifflin  [24] 
and  Stoer  [28]  used  the  active  set  strategy  in  their  algorithms  to  solve  a constrained 
least-squares  problems.  Fletcher  and  Jackson  [8],  Gill  and  Murray  [9]  designed  their 
algorithms  based  on  the  active  set  strategy  to  solve  a more  general  quadratic  pro- 
gramming problem  with  symmetric  Hessian  matrix  and  inequality  constraints.  More 
and  Toraldo  [25]  combine  active  set  strategy  with  the  gradient  projection  method  of 
Goldstein  [11]  and  Levitin  and  Polak  [22]  to  solve  a quadratic  programming  problem 
with  box  constraints. 

1.3  Dual  Active  Set  Algorithm 

William  W.  Hager  presented  a new  active  set  algorithm,  the  Dual  Active  Set 
Algorithm  [15],  for  solving  a minimization  problem  with  equality  constraints  and  box 
constraints.  The  dual  active  set  algorithm  maximizes  the  corresponding  Lagrangian 
dual  problem  with  an  unconstrained  Lagrangian  multiplier  vector.  Starting  from  an 
arbitrary  Lagrangian  multiplier,  it  identifies  a set  of  active  bound  constraints  which 
is  called  an  active  index  set.  The  dual  active  set  algorithm  obtains  a solution  to 
the  corresponding  Lagrangian  dual  problem  in  a finite  number  of  iterations.  The 
procedure  for  changing  the  active  index  set  in  this  algorithm  is  the  inverse  of  the 
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procedure  of  the  primal  active  set  strategy.  In  this  algorithm,  indices  are  dropped 
out  from  the  active  index  set  while  maximizing  the  dual  problem. 

We  briefly  describe  the  Dual  Active  Set  Algorithm  as  follow: 

Consider  the  following  minimization  problem, 

min  /(x) 

s.t.  h(x)  = 0 (1.2) 

l < x < u 

where  / : RJ1  — ► R and  h : Rn  — ► Rm  are  continuously  differentiable  functions,  / and 
u are  constant  vectors. 

The  unconstrained  Lagrangian  dual  problem  of  the  minimization  problem  (1.2) 
with  the  Lagrangian  multiplier  A 6 Rm  is  obtained  by  relaxing  the  equality  constraint 
h(x)  = 0 and  is  defined  as 

max  L( A)  s.t.  A G Rm  (1-3) 

where 

L(X)  = inf  L(\,x)  s.t.  I < x < u (1.4) 

and 

L(  A,  x)  = f(x)  + A Th(x)  (1.5) 

Some  notation  and  auxiliary  problems  should  be  introduced  before  we  introduce  the 
algorithm. 

• Let  1b  be  a subset  of  the  index  set  {1,2,  ...,n}. 

• Let  xb  denote  the  vector  with  components  x,  associated  with  indices  * € Ib- 


Lb(X)  = inf  L( A,  x)  s.t.  lB  < xB  < uB 


(1.6) 
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• Given  a vector  z € -ft”,  define 

LZB(\)  = inf  L(\,x)  s.t.  xB  = zB  (1.7) 

The  following  algorithm  is  given  in  [15].  It  describes  the  procedure  to  compute  the 
next  iteration  if  A*  is  the  current  iteration. 


Dual  Active  Set  Algorithm 


Given  A*,  let  j = 0,  i/0  = A*,  and  define 


where 


I B0  = {i:  Zi  = li  or  z>  = u,} 


z = x(Afc)  = arg  min  L( A*,  x)  s.t.  I < x < u. 


Sub-iteration:  Let  fij  maximize  (A)  over  A and  define  g,(t)  = vj  + t(fij  — i/j). 
Determine  the  largest  interval  [0,  t\,  t > 0,such  that 
£&,(/*(*))  = Lbi(p(t))  for  every  t € [0,/]. 

If  t < 1,  then  put  i/j+i  = g,(t).  The  set  7bj+1  is  obtained  by  deleting 
from  IB}  those  indices  i £ IB.  with  the  property  that 

L , =n 

dxi  IA=A,x=x 

where  x is  a minimizer  in  ( 1.7)  associated  with  A = A = g,(t). 
Increment  j and  repeat  the  sub-iteration. 

If  t > 1,  then  set  A*+1  = fij,  increment  k,  and  proceed  to  the  next 


iteration. 


8 


The  dual  active  set  algorithm  identifies  a set  of  active  bound  variables  while 
minimizing  the  problem 

L(Xk , x ) s.t.  I < x < u 

with  a given  multiplier  A*.  The  constraint  h(x)  = 0 is  relaxed  until  termination 
while  a set  of  bound  constraints  is  maintained  and  constantly  changed  throughout 
the  iterations.  In  the  successive  iterations  which  maximize  LB( A),  some  of  the  active 
bound  constraints  will  be  dropped  to  increase  the  dual  objective  function  value.  Since 
the  algorithm  constantly  increases  the  dual  function  value  while  maintaining  a set  of 
bound  constraints,  the  set  of  active  bound  constraints  will  never  repeat.  With  only 
finitely  many  choices  for  the  different  active  index  sets,  the  dual  active  set  algorithm 
will  attain  a solution  of  the  Lagrangian  dual  problem  (1.3)  in  finitely  many  iterations 
which,  at  the  same  time,  also  solve  the  minimization  problem  (1.2). 

The  following  theorem,  which  shows  the  finite  termination  and  convergence  prop- 
erty of  the  dual  active  set  algorithm,  can  be  found  in  [15]. 

Theorem  1.3.1  If  f and  h are  continuously  differentiable  on  Rn , L satisfies  the  strong 
convexity  assumption 

L{\y)  > L(\,x)  + yxL(\,x)(y  - x)  + a(y  - xf 

for  a constant  a > 0,  and  there  exists  a maximizer  pj  of  LBj  for  each  j , then  the 
dual  active  set  algorithm  reaches  a solution  of  (1.3)  in  a finite  number  of  iterations 
and  sub-iterations. 

A computational  study  of  large-scale  quadratic  network  problems  comparing  var- 
ious algorithms  to  the  dual  active  set  algorithm  can  be  found  in  Hager  and  Hearn 
[16].  Their  numerical  experiments  have  shown  favorable  results  for  this  method.  It 
is  also  pointed  out  in  that  paper  that  the  dual  active  set  algorithm  is  very  effective 
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for  the  quadratic  network  problems  if  the  test  problems  are  preprocessed  by  a few 
number  of  conjugate  gradient  iterations. 


CHAPTER  2 

DUAL  ACTIVE  SET  ALGORITHM  FOR  LINEAR  PROGRAM 


2.1  Introduction  and  Review 

In  this  chapter,  we  consider  the  following  linear  programming  problem  (LP) 

min  cTx 

s.t.  Ax  — b (2.1) 

l < x < u 

where  c,  /,  u and  x are  vectors  in  Rn,  b is  a vector  in  Rm,  and  A is  a m x n matrix. 

Linear  programming  has  been  used  extensively  in  a variety  of  areas  in  applied 
mathematics,  business,  industry,  military,  etc.  Among  those  methods  developed 
in  the  past  half  century  for  solving  the  linear  programming  problems,  the  Simplex 
method  and  various  types  of  interior  point  methods  are  most  noticeable  from  both 
the  theoretical  and  computational  points  of  view.  The  Simplex  method,  developed 
by  George  B.  Dantzig  in  1947,  starts  from  a vertex  of  the  polyhedron  and  then  travels 
along  the  edges  from  one  vertex  to  another  to  minimize  the  problem.  This  method  had 
been  used  in  linear  programming  without  serious  competition  until  the  introduction 
of  the  idea  of  the  interior  point  method  applied  to  linear  programming  by  Narendra 
Karmarkar  in  1984,  see  Bazaraa,  Jarvis,  and  Sherali  [2].  For  an  introduction  to 
various  types  of  interior  point  methods,  see  Wright  [30].  Many  appealing  proofs  of 
termination  in  polynomial  time  for  various  types  of  interior  point  methods  have  been 
extensively  studied  in  the  past  decade.  A collection  of  journal  articles  and  research 
reports  on  the  interior  point  methods  can  be  seen  in  Kranich  [21]. 

Ever  since  the  development  of  the  interior  point  methods  for  linear  programming, 
numerical  experiments  on  comparison  of  the  Simplex  method  and  interior  point  meth- 
ods have  been  done  by  many  people.  For  examples,  see  Adler  et  al.  [1],  Tomlin  [29]. 
Results  so  far  show  that  interior  point  types  methods  behave  favorably  on  large  linear 
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programming  problems.  However,  in  spite  of  the  presence  of  a worst-case  complexity 
example  for  the  Simplex  method  which  travels  through  every  vertex  of  the  constraint 
polyhedron,  see  Klee  and  Minty  [20],  the  Simplex  method  still  shows  its  strength  on 
solving  smaller  problems. 

The  general  linear  programming  problem  (LP)  can  also  be  solved  by  making 
suitable  modifications  to  the  dual  active  set  algorithm  presented  in  section  (1.3). 
This  is  the  objective  of  this  chapter. 

In  the  remainder  of  this  chapter,  we  introduce  the  dual  active  set  algorithm  for 
the  linear  programming  problem.  Modifications  of  the  dual  active  set  algorithm  are 
described  in  detail.  A variation  of  this  algorithm  in  the  procedure  of  adding  active 
bound  constraints  is  given  in  section  2.2.  This  method  is  sound  provided  a certain 
assumption  is  imposed  on  this  (LP)  problem.  Convergence  properties  are  studied 
through  a series  of  lemmas  and  theorem  in  section  2.3.  Numerical  implementation  of 
the  test  problems  from  Netlib  is  discussed  in  section  2.4.  Also,  numerical  experiment 
results  axe  presented  in  the  end  of  this  chapter. 

2.2  The  Algorithm  and  its  Variation 

The  goal  of  this  section  is  to  present  the  dual  active  set  algorithm  for  the 
linear  programming  problem  (LP).  This  algorithm  is  derived  from  the  dual  active  set 
algorithm  in  section  (1.3).  We  present  two  versions  for  this  algorithm.  The  first  one 
does  not  need  any  prior  assumption  to  the  linear  programming  (LP)  while  the  second 
one  does  need  one. 

Consider  the  Lagrangian  L( A,  x ) of  (LP)  with  a multiplier  vector  A 6 R171 

L(X,x)  = cTx  + \T(b-Ax) 

= A Tb  + (c  — At\)tx 
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It  is  well  known  that  the  linear  programming  problem  (LP)  is  convex  but  that  it 
does  not  satisfy  the  strong  convexity  assumption  demanded  in  the  Dual  Active  Set 
Algorithm  for  the  general  problem  described  in  the  previous  chapter.  That  is,  there 
exists  a constant  a > 0 such  that 

L(X,  V ) > L(X,  x)  + S7xL(X,  x)(y  - x)  + a(y  - x)2, 

for  all  x and  y in  the  feasible  domain.  In  fact,  the  Lagrangian  L( A,  x)  for  (LP)  only 
satisfies 

■^(A,  J/)  — L(X,x)  T V*L(A,  x)(y  x), 

for  all  x and  y in  the  feasible  domain.  Under  the  framework  of  the  Dual  Active 
Set  Algorithm,  we  modify  this  algorithm  to  overcome  the  problem  caused  by  the 
lack  of  the  strong  convexity  assumption  when  we  apply  this  algorithm  to  the  linear 
programming  problem.  In  this  section,  modification  of  the  original  Dual  Active  Set 
Algorithm  is  analyzed  and  discussed.  For  consistency,  we  use  the  same  notations  for 
the  Dual  Active  Set  Algorithm  as  we  do  in  the  linear  programming  problem. 

The  following  is  a list  of  expressions  for  the  notations  and  auxiliary  problems  used 
in  the  dual  active  set  algorithm  in  section  (1.3)  corresponding  to  that  of  the  linear 
programming  problem  (LP): 

1.  Define  /#,  In,  x#,  and  x n the  same  as  in  section  (1.3). 

2.  Denote  A,  as  the  *-th  column  of  the  matrix  A. 

3.  Given  an  index  set  Ib  and  a corresponding  index  set  In,  denote  the  matrix  B 
as  a sub-matrix  of  the  matrix  A with  columns  of  A associated  with  indices  in 
Ib-  Similarly  for  N. 
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4.  For  A £ Rm , define 


L{ A)  = inf  L( A,x) 

s.t.  1 < x < u 

= XTb  + inf  (c  — AT X)T x 

s.t.  1 < x < u 

5.  Given  an  index  set  Ib , define 

Lb( A)  = inf  L( A,  x) 

s.t.  Ib  < xb  < ub 

= XTb  + inf  (c  — ATX)Tx 

s.t.  Ib  < xb  < ub 

6.  Given  a vector  z,  and  an  index  set  7s,  define 

Lzb{ A)  = inf  L(X,x) 

s.t.  xB  = zB • 

= XTb  + inf  (c  — AT X)Tx  s.t.  xb  = zb- 

Note  that  is  unconstrained  in  Lb( A)  and  LZB(X). 

The  Dual  Active  Set  algorithm  is  designed  to  solve  the  Lagrangian  dual  problem 
of  (LP) 

max  L{ A) , A £ Rm. 


The  dual  approach  of  this  algorithm  has  an  advantage  that  the  Lagrangian  multiplier 
is  not  restricted  while  the  disadvantage  is  the  lack  of  smoothness  in  the  objective 
function.  In  the  following,  we  expose  the  problems  that  the  original  dual  active  set 
method  will  encounter  through  one  iteration  if  we  employ  it  on  (LP)  directly. 

At  the  beginning  of  the  A;-th  iteration,  A*  is  known,  we  solve  the  minimization 
part  of  the  dual  problem.  A solution  to 

min  ( c — ATXk)Tx 


s.t.  I < x < u 


(2.2) 


can  be  determined  by  the  following  simple  rules: 


1.  If  (c  — AT\k)i  > 0 then  x,  = 
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2.  If  (c  — ATXk)i  = 0 then  / < x < u. 

3.  If  (c  — Ar\k)i  < 0 then  x,  = u,-. 

therefore  the  active  index  set  Ib  and  the  non-active  index  set  In  could  be  set  as 
Ib  = {i  : [c  - ATA*], 0 },  In  = {l,2,...,n}  - IB , 

as  well  as 

zb  = {zi  ■ Zi  £ z,  i E Ib  , z = x(Afc)  = arg  min  L(\k,x)  s.t.  I < x < u}. 
Since  z solves 

min  L(Xk,  x)  s.t.  I < x < u, 
z also  solves  the  following  two  problems: 

1.  min  L(Xk,  x)  s.t.  Ib  < xb  < 

2.  min  L(Xk,x)  s.t.  xb  = zb 

therefore,  we  have  Ls(Xk)  = LzB{Xk)  = L(Xk).  This  makes  the  value  of  dual  function 
equal  to  that  of  those  two  auxiliary  problems  at  the  beginning  of  each  main  iteration. 

In  the  algorithm  we  maximize  LB(X)  through  a sequence  of  sub-iterations  by 
finding  proper  directional  vector  and  corresponding  step  size  along  the  direction. 

If  (c  — ATX)i  0,  for  some  i € In,  then 

Lzb( A)  = XTb+ inf(c  — ATX)Tx  s.t.xs  — ZB 

(2.3) 

= — oo 

because  xn  is  unconstrained.  Therefore  we  cannot  use  the  gradient  of  LB(Xk)  directly 
as  the  directional  vector  to  maximize  the  problem.  Also,  caution  should  be  taken  to 
choose  a proper  step  size  along  the  directional  vector  which  will  change  the  current 
active  index  set. 
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Now  we  discuss  how  to  find  a suitable  direction  vector  and  a step  size  along  this 
vector  such  that  we  can  maximize  the  dual  problem  and  also  maintain  a set  of  active 
indices.  Moreover,  the  step  size  chosen  may  drop  out  some  indices  from  the  current 
active  index  set.  Let  d be  the  directional  vector  and  a be  the  corresponding  step 
size.  The  strategies  we  used  to  overcome  these  problems  are 

1.  Find  d such  that  along  this  directional  vector  it  increases  the  value  of  LB  and  also 

maintains  the  current  active  index  set. 

2.  Find  step  size  a such  that  LzB(\k  + sd)  = LB{\k  + sd ) for  0 < s < a. 

To  find  a suitable  directional  vector  d that  satisfies  the  above  statements,  we 
consider  the  following  sub-iterations. 

Let  j = 0,  Po  = A*,  and  /3(a)  = (3j  + adj , for  a > 0.  where  dj  is  determined  to 
satisfy  the  following  two  conditions: 

1\  /3(a)  maximize  LB(/3(a))  and  (c  - v4T(/?(a)));v  = 0. 

2’.  LB((3(a))  = Lb((3(o))  for  q < a where  a is  a positive  number. 

Note  that  LB(/3(a))  and  Lb(/3(o))  can  be  written  as 

^b(^(q))  = Pjb  + ad%b  + inf[(c  - AT/3j  - aATdj)Tx]  s.t.  xB  = zB. 

Lb((3( a))  = f3jb  + ad/jb  + inf[(c  - AT /3j  — aATdj)Tx\  s.t.  lB  < xb  < Ub- 

Since  these  two  auxiliary  problems  differ  only  on  the  minimization  part,  the  di- 
rectional vector  dj  can  be  chosen  in  a way  to  satisfy  both 

1 ”■  [dj A]i  = 0,  for  i € In,  that  is  dj N = 0. 
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2”. 


dj (b  — Ax)  > 0 
=>•  dj (b  — Bxb  — Nxn)  > 0 

=>•  djb  — djBxB  — dj Nxn  > 0 since  dj N = 0 
=>  dj (b  — Bxb)  > 0 

A directional  vector  satisfying  the  above  requirements  may  be  chosen  as  the  projection 
of  the  gradient  ( b — Bzb)  of  LB  onto  the  null  space  of  NT . 

After  we  find  the  directional  vector  dj,  a line  search  should  be  performed  to  obtain 
a new  Lagrangian  multiplier.  Exact  line  search  is  used  in  this  algorithm.  A step  size 
is  determined  by  finding  the  smallest  step  such  that  one  or  more  components  of  the 
coefficient  vector  ( c—ATf3j—aATdj ) of  the  minimization  part  of  the  auxiliary  problem 
Lzb  equal  zero.  This  exact  line  search  enables  us  to  drop  some  indices  from  the  set 
Ib  to  the  set  1^. 

Proceeding  to  the  sub-iteration,  this  algorithm  will  encounter  two  different  pos- 
sible cases  with  respect  to  the  directional  vector  and  the  corresponding  step  size. 
During  the  sub-iterations,  the  directional  vector  d could  have  two  cases: 

I.  If  d ^ 0,  then  proceed  to  find  the  step  size. 

(a) .  IfO  < a < oo,  active  indices  may  be  dropped  from  Ib  to  /at. 

(b) .  If  a = oo,  we  can  show  that  this  problem  is  infeasible. 

II.  If  d = 0,  then  we  cannot  increase  LB  under  the  current  active  set  anymore,  A 

new  direction  can  be  found  by  solving  the  following  least-squares  problem 

min||Ar  — 6||  s.t.  xb  = zb,  l < x < u. 
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Let  y be  a solution  of  the  above  least  square  problem. 

If  Ay  = b , then  we  can  show  that  y is  a solution  to  (LP). 

Otherwise,  we  add  in  those  indices  which  have  yt  = /,  or  y,  = tx,,  for  i E In  to 
the  active  index  set.  Then,  we  search  along  the  direction  ( b—Ay ) and  determine 
a proper  step  size  such  that  some  indices  can  be  dropped  out  from  the  current 
active  index  set.  This  finishes  a main  iteration  and  we  can  show  that  the  dual 
objective  has  been  strictly  increased. 

Now  we  present  the  dual  active  set  algorithm  for  the  linear  programming  problem. 
Note  that  the  sub-iterations  between  the  main  iteration  A*  and  A k+i  are  denoted  by 

ft- 


Dual  Active  Set  Algorithm  for  Linear  Programming  (I) 


Main  iteration: 
Given  A*, 

Let 


2 € arg  min  {crx  + Af  (6  - Ax)  : l < x <u}, 


and 


Ib  = {*’  : (c-  AT Xk)i  # 0}. 


Sub-iteration: 


Let  j = 0,  /?o  = A*,  IBo  = IB,  (3(a)  = (3j  + adj. 

dj  : the  projection  of  (6  — BjZBj)  into  the  null  space  of  Nj. 
If  dj  ^ 0,  then 

find  aj  such  that  LB]((3(a))  = LzB](/3(a)),  a < aj,  and 
[c  — AT/3(aj)\i  = 0 for  some  i 6 IB]. 


(2.4) 
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If  oij  = oo,  then  (LP)  is  unbounded,  stop, 
otherwise, 

jb]+1  = Ib}  ~ {i  € IB],  [c  — AT (3 (a j)\i  = 0}. 

Pi+i  = /%;)> 

proceed  to  the  next  sub-iteration, 
otherwise 

Let  A k be  the  last  f3j , /g . be  the  last  active  index  set 

and  let  y denote  a solution  to  the  least-squares  problem  (LS), 

min  ||Ar  — 6||  s.t.  xg.  = zg.,  / < x < u (LS) 

If  Ay  = 6,  then  y is  a solution  to  (LP).  Stop, 
otherwise, 

update  7g.  according  to  y,  define  A(a)  = Ajt  + a(b  — Ay),  let  a be  the 
smallest  a > 0 such  that  [c  — ATA(a)],  = 0 for  some  i E Ib *• 

If  a = oo,  then  (LP)  is  unbounded,  stop, 
otherwise, 

and  Ajfc+i  = A(a),  proceed  to  next  main  iteration. 


We  now  consider  another  version  of  this  algorithm  when  an  additional  assumption 
is  added  to  the  (LP)  problem.  The  assumption  is: 

Assume  that  the  columns  of  the  non-active  sub-matrix  N of  the  constraint  matrix  A 
are  linearly  independent  in  each  sub-iteration. 

A formula  for  each  directional  vector  d can  be  set  up  in  the  following  way: 

Since  d is  chosen  as  the  projection  of  (6  — Bzb)  onto  the  null  space  of  NT , for  some 
vector  v,  we  can  derive  a formula  for  d as  follows:  Write  the  vector  (b  — Bzb)  as  a 
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direct  sum  of  d and  Nv, 

( b - Bzb)  = d + Nv.  (2.5) 

Multiplying  both  sides  of  the  equation  2.5  by  NT,  we  get 

NT(b-BzB ) = NTd  + NTNv 

=>•  NT(b  — Bzs ) = NtNv  (sinceNTd  = 0.)  (2-6) 

By  the  assumption  that  the  columns  of  N are  linearly  independent,  it  can  be  shown 
that  (ArTAr)_1  exists.  Hence,  by  multiplying  (jVr./V)_1  to  both  sides  of  equation  2.6, 
the  vector  v can  be  expressed  as 

u = (NT N)'1  NT (b  - Bzb ) 

By  substituting  the  vector  v back  into  the  equation  (2.5),  we  obtain  a formula  for 
the  directional  vector  d,  namely, 

d = (/  - Ni^Ny'N^ib  - Bzb)  (2.7) 

If  the  directional  vector  dj  is  nonzero,  we  can  show  that  this  dj  is  an  ascend  direction 
for  LzB(Xk)  with  a given  A k,  otherwise,  if  dj  = 0,  then  by  the  projection  equation,  we 
have  Nv  = (6  — Bzb)  for  some  vector  v.  We  discuss  two  cases  which  lead  either  to 
an  optimal  solution  for  this  dual  problem  or  to  an  update  of  the  current  active  index 
set  and  then  proceed  to  the  next  main  iteration. 

I.  If  d ^ 0,  use  the  same  process  in  algorithm  (I). 

II.  If  d = 0,  by  the  projection  formula  ( b — Bzb)  = d 4-  Nv^,  we  have  an  equality 

equation  Nv n = (b  — Bzb ) with  an  unknown  vector  ujy.  The  least-squares 
solution  to  this  equation  is 


vN  = {NTN)~1NT{b-BzB). 


(2.8) 
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(%>• 

(2). 


If  bv  < vn  < u^v,  we  can  show  that  [zg,  vfj]T  is  a minimizer  for  (LP). 


If  Vi  > Ui  or  Vi  < /,•  for  some  i € In,  then  an  update  on  the  active  index 
set  is  required.  We  add  in  those  indices  IB\  = {*  : u,-  > u,-  or  u,  < /,}  to 
the  current  active  index  set  and  also  update  z by 

Ui  ; if  Vi  > Ui 


Zi  = 


i G /#!. 


, j if  ^ I« 


Set 


d = [/  - W(iVrAT)-1WT](6  - Bzb) 


Delete  those  indices  i G IB 1 satisfy  one  of  the  following 
(a),  dfbi  < 0 if  Zi  = u,-. 

<Fbi  > 0 if  z,-  = 

where  6^  is  a column  of  B1.  Update  IB,  In,  and  zB,  repeat  until  no  indices 
in  IB i satisfy  the  conditions  (a)  and  ( b ). 

This  process  will  be  terminated  in  finitely  many  iteration  since  IBi  is  a 
finite  set.  We  can  show  this  procedure  increases  the  dual  function  value 
and  also  updates  the  active  index  set. 


We  denote  this  version  as  algorithm  (II)  and  present  it  as  follow: 


Dual  Active  Set  Algorithm  for  Linear  Programming  (II) 

Given  A0, 
let 

IB  = {i  : [c  — AT\k]i  7^  0 },  IN  = {l,2,...,n}  - IB, 

where 

zb  = {zi  : Zi  € z,i  € IB  where  z = x(Ajt)  = argminL(\k , x ) s.t.  I < x < u.} 
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Sub-iteration: 

Let  j = 0,  Po  = \k,  Find  dj  such  that 

dj  = [I-  JV^TV)"1  jVT](& - Bzb) 


If  dj  ^ 0, 

Let  (3(a)  = Pj+adj.  Find  the  largest  step  size  aj  such  that  Lb(P(q))  = LzB(P(a )) 
for  all  a < aj  and  [c  — ATP(aj)]i  — 0,  for  some  i € IB-  Delete  i from  IB,  update  1^ 
and  zB ■ Set  /?J+1  — pj  -f  ajdj,  proceed  to  the  next  sub-iteration. 

If  dj  = 0, 

Solve  the  linear  system  Nv  = b — Bzb  for  v. 


1.  If  / < v < u,  stop,  x*  = [zg,  vT]T  is  a minimizer  for  (LP). 

2.  If  Vi  > Ui  or  Vi  < li,  add  these  indices  IB i = [i  : u,-  > u,  or  u,-  < /,}  to  the 
current  active  set,  update  zB  by  setting 

Ui  ; i f Vi  > Ui 


Set 


if  Vi  < f 


i € Ib1- 


d = [I  — N(NTN)-1NT](b  - Bzb) 


Delete  those  indices  i € IB i and  satisfy  one  of  the  following 


(a)  cFbi  < 0 if  Zi  = Ui. 

(b)  dfbi  > 0 if  Zi  = 


where  fej  is  a column  of  B1.  Update  IB , /jv>  and  zB,  repeat  until  no  indices  in 
IB 1 satisfy  the  conditions  (a)  and  (b). 

Let 

= Pj  + oid 
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denote  this  new  active  set  as  Ig,  update  In  and  zg,  proceed  to  the  next  itera- 
tion. 


At  each  sub-iteration,  a column  from  the  active  sub-matrix  B may  be  added 
into  the  non-active  sub-matrix  A.  Therefore,  A is  constantly  extended  until  d = 
0.  However  the  directional  vector  d = (I  — A(ArA)-1  NT)(b  — Bzg)  need  not  be 
recomputed  entirely  whenever  A is  changed.  Since  we  only  add  in  one  column  to  A 
at  each  sub-iteration,  the  directional  vector  d can  be  updated  efficiently  by  applying 
a rank  one  correction  procedure  to  the  inverse  of  the  old  ( NTN ).  To  find  the  new 
(ArA)_1,  which  needs  the  most  computational  efforts  in  updating  d,  an  exercise  in 
Luenberger  [23]  page  361,  shows  that 


Proposition  2.2.1  Let  A = [a,  A]  be  a matrix  constructed  by  adding  a column  vector 
a to  a matrix  N,  then  (ArA)_1  can  be  found  from  (ATA)-1  by  the  formula 

e -taTN(NTN)~1 

—e(NTN)~1NTa  (ArA)-1[/  + eNTaaTN(NT A)"1]  _ 

where 


(ArA)-1  = 


t = 


aTa  - arA(ATA)_1ATa 


and 


AJA  = 


cr  a 


aTN 


NTa  NtN 


Proof  of  Proposition  2.2.1  If  the  inverse  matrix  of  (ArA)  exists,  it  can  be  written 
as 

where  0 ^ a £ R,  (I  is  a column  vector,  and  C is  a square  matrix. 
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Since  ( NTN ) • ( NTN ) 1 = /,  we  have 

aTaa  + aTN/3  = 1 

(2.9) 

aTaPT  + aT NC  = 0 

(2.10) 

NTa/3T  + NtNC  = I 

(2.11) 

Multiplying  (2.11)  by  ( NTN ) 1 gives 

c = (NTN)~l  - (Nt N)-1  Nt a(3T . 

(2.12) 

Substituting  (2.12)  into  (2.10),  we  get 

T _ -a't'NiNT'N)-1 
**  ~ aTa-aTN(NTN)-1NTa 

Substituting  (3T  into  (2.9)  and  (2.12)  respectively  we  have 

1 

a~  aTa-  aI’^(ATTAf)-17VTa 

and 

C = (A^r7V)_1[/  + aNTaaTN(NT  N)'1]. 

□ 


2.3  Convergence 

We  establish  the  convergence  of  the  dual  active  set  algorithm  for  linear  pro- 
gramming in  this  section.  Convergence  properties  for  both  versions  will  be  discussed 
respectively.  The  goal  is  to  show  that  the  dual  active  set  algorithm  for  linear  pro- 
gramming attains  an  optimal  solution  to  (LP)  in  finitely  many  iterations.  We  reach 
this  goal  through  a sequence  of  lemmas.  Note  that  lemma  2.3.2,  2.3.3,  2.3.4,  and 
theorem  2.3.1  can  be  found  in  Hager  and  Shih  [18]. 
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Lemma  2.3.1 

(1) .  If  the  directional  vectors  defined  by  the  projection  of  the  gradient  of  LB(\k)  with 

respect  to  x is  nonzero,  then  it  is  an  ascend  direction  for  the  dual  problem. 

(2) .  If  Ay  ^ b where  y is  a solution  to  the  least-squares  problem  (LS),  then  the 

directional  vector  (b  — Ay)  used  in  adding  new  indices  to  the  active  index  set  is 
an  ascend  direction  for  the  dual  problem. 

Proof  of  Lemma  2.3.1 

( 1).  Since  the  directional  vector  dj  is  obtained  by  projecting  the  gradient  ( b—BjZBj ) 
of  LzB}( A*)  onto  the  null  space  of  NT , the  vector  [(b—BjZBj)—dj]  is  perpendicular 
to  dj.  We  have 

dJ(b-BjZBj)  = <fj[{b  - BjZBj)  - d:  + dj]  (2.14) 

= Mil2- 

Since  z is  a solution  to 


min  cTx  + A f (6  — Ax)  s.t.  I < x < 


u , 


we  have 


(c  — At fij)  > 0 if  Zi  — /, 

(c  - AT ft j)  < 0 if  Zi  = Ui. 

Hence,  by  the  above  relations  and  (NTdj)ff  = 0,  we  have 

L(P(a)) 

= L(()j  + adj) 

= fijb  + adj b + min(c  — AT fij  — aATdj)Tx , s.t.  I < x < u 
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= (3jb  + adj b+  (c  — Atj3j)zb  - a{ATdj)Tz 
= f3jb  + adj b — a(BTdj)T zb  + min(c  — AT/3j)Tx,  s.t.  I < x <u 
= L((3j)  + adj b - ad f Bzb 
= L((3j)  + adj ( b - Bzb ) 

= L(f3j ) + a\dj\2,  for  a > 0. 

If  dj  0,  then  <f-(b  — Bzb ) = \dj\ 2 > 0.  Therefore,  with  a > 0,  the  above 
relation  shows  that  L(fi(a))  > This  proves  that  d is  an  ascend  directional 

vector. 

(2).  If  Ay  ^ b and  y is  a solution  to  the  least-squares  problem  (LS),  then  , by  the 
first-order  necessary  condition  for  the  solution  y of  (LS),  we  have  the  following 
relations: 


{Ay  - b)TAi  = 0 if  * € JV 

and  /,•  < j/i  < u, 

(2.15) 

{Ay  - b)T Ai  >0  if  i G N 

and  y,  = /,• 

(2.16) 

{Ay  — b)T  Ai  < 0 Hie  N 

and  y,  = u, 

(2.17) 

By  the  definition  of  the  active  index  set  B, 

(c  - ATXk)i  = 0, 

i e N. 

(2.18) 

Hence,  by  the  relations  (2.15)  ~ (2.18),  we  got 

(c  - ATA(a))i  = 0 if  ieN 

and  /,  < yi  < Ui 

(2.19) 

(c  - ArA(a)),-  >0  if  ieN 

and  y,  = /, 

(2.20) 

(c  - >lTA(a)),'  <0  if  * 6 N 

and  y,  = u, 

(2.21) 

for  a > 0. 

However,  since  ys  = zb  and  z is  a solution  to 


min  cTx  + A^(6  - Ax)  s.t . I < x < u,  (2.22) 
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we  have 

(c  — i4rA(a)),-  >0  if  i € B and  y,  = z,-  = L (2.23) 

(c  — ATA(a)),  <0  if  i € B and  y,-  = = iti  (2.24) 

for  a > 0. 

Therefore, 

L(X(a)) 

= L(Xj  + a(b  - Ay)) 

= A J b + q(6  - Ay)Tb  + min(c  - AT\k  - aAT(b  - Ay))T  x s.t.  I < x <u 

= \kT  b + (c  — AT\k)T  z + a(b  — Ay)Tb  — a(b  — Ay)T Ay 

= \kTb+(c-  AT \k)T z + a\\b  - Ay\\2 
= L(Xk)  + o||6 - Ay\\2 

Since  the  step  size  a > 0,  and  Ay  ^ b implies  ||6  — Ay\\2  > 0,  these  show  that 
( b — Ay)  is  an  ascend  direction  for  the  dual  problem. 

□ 

Lemma  2.3.2  If  the  step  size  a is  infinite , then  the  linear  programming  problem  (LP) 
is  infeasible. 

Proof  of  Lemma  2.3.2  There  are  two  places  where  a is  involved,  within  the  sub- 
iteration  and  after  solving  the  least  squares  problem  (LS)  which  is  (1.)  and  (2.)  in 
lemma  2.3.1. 

(1.)  Within  the  sub-iteration:  From  lemma  2.3.1  we  have 


L(fi{a))  = L(fij)  + a\di \2 


(2.25) 
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(2.)  After  solving  the  least-squares  problem  (LS):  From  lemma  2.3.1  we  have 

L( A(a))  = L(Xk)  + a\\b  - Ay\\ 2 (2.26) 

If  the  step  size  a,  in  either  equation  (2.25)  or  (2.26),  becomes  unboundedly  large 
then  the  dual  objective  function  value  tends  to  infinity.  By  the  Weak  Duality  Theory, 
the  optimal  objective  value  of  the  primal  problem  is  greater  than  or  equal  to  the 
optimal  objective  value  of  the  dual  problem.  Therefore,  the  primal  problem  (LP)  is 
infeasible.  □ 

Next  we  examine  the  optimal  condition  for  termination  of  this  algorithm. 

Lemma  2.3.3  If  Ay  = b where  y is  a solution  to  the  least  squares  problem,  then  y is 
an  optimal  solution  to  the  linear  programming  problem  (LP)  and  Xk  is  optimal  in  the 
dual  problem. 

Proof  of  Lemma  2.3.3  If  y is  a solution  to  (LS)  and  Ay  = b,  then  y is  feasible  for  the 
primal  problem  (LP).  Since  A*  is  unconstrained,  it  suffices  to  show  that  the  value  of 
the  primal  equals  the  value  of  the  dual. 

Since  z is  chosen  to  minimize  the  problem  (2.22),  we  have 

(c  - ATXk)i  = 0 if  /,  < y,  < Ui 
(c  - ATXk)i  > 0 if  yi  = U 
(c  - ATXk)i  < 0 if  yi  = u, 

and  ys*  = 2b*.  we  have 

L(Xk))  = cTz  + Xl(b-Az ) 

= (c-  ATXk)Tz  + Xlb 
= (c-  ATXk)g.zB*  + Xjb 
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= (c  - AT\k)Ty  + XTkb 
= cTy  + \l(b-  Ay) 

= cTy 

where  IB*  is  the  last  active  index  set.  This  shows  that  y is  optimal  in  the  primal  and 
A k is  optimal  in  the  dual.  □ 

Lemma  2.3.1  Let  IBi  denote  the  last  active  index  set  ati-th  main  iteration.  IB<  ^ IB3 

ifi 

Proof  of  Lemma  2.3.4  Since  the  step  size  chosen  at  each  iteration  is  strictly  posi- 
tive, equation  (2.25)  and  (2.26)  imply  that  the  dual  objective  function  value  at  each 
iteration  increases  strictly.  Therefore,  the  active  index  set  LB,  never  repeats.  □ 

Theorem  2.3.1  In  a finite  number  of  iterations,  the  dual  active  set  algorithm  for  lin- 
ear programming  either  determines  that  (LP)  is  infeasible,  or  it  obtains  an  optimal 
solution. 

Proof  of  Theorem  2.3.1  If  the  linear  program  (LP)  has  a feasible  solution,  then  the 
step  sizes  a that  appear  in  the  dual  active  set  algorithm  are  finite  and  strictly  positive. 
Thus,  the  dual  objective  function  value  at  each  iteration  is  strictly  increased.  Since 
each  iteration  corresponds  to  a different  active  index  set  by  lemma  2.3.4,  and  there  are 
only  finitely  many  combinations  of  active  index  sets.  The  dual  active  set  algorithm 
for  linear  programming  will  attain  an  optimal  solution  of  (LP)  in  a finite  number  of 
iterations.  □ 

With  the  assumption  added  in,  we  examine  the  convergence  property  for  algorithm 
(II).  Without  loss  of  generality,  we  also  assume  (LP)  has  a feasible  solution. 
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Theorem  2.3.2  Assume  that  the  columns  of  N are  linearly  independent  at  each  it- 
eration of  the  dual  active  set  algorithm  for  (LP),  and  (LP)  has  a feasible  solution. 
The  algorithm  (II)  reaches  a solution  of  (LP)  in  a finite  number  of  iterations  and 
sub-iterations. 

Proof  of  Theorem  2.3.2  With  the  assumption  that  the  columns  of  N are  linearly 
independent  at  each  sub-iteration,  we  first  show  that  the  number  of  sub-iterations 
in  each  main  iteration  should  be  less  than  or  equal  to  m.  Since  each  sub-iteration 
adds  a new  column  to  the  matrix  N,  the  matrix  N is  constantly  extended.  When 
the  number  of  columns  of  N is  m,  N is  a m x m non-singular  matrix.  This  implies 
I = N(NT N)*1  NT , hence  d = 0.  Therefore,  the  number  of  sub-iterations  in  each 
main  iteration  must  be  less  than  or  equal  to  m. 

If  dj  = 0 and  / < v < u,  x*  = [zB,  vT]T  satisfies  the  constraints  in  (LP).  Moreover, 
it  satisfies  the  optimality  condition  that 

1.  [c  — AT/3j\i  < 0 if  Xi  = Ui. 

2.  [c  — AT/3j]i  > 0 if  Xi  = L. 
therefore  x*  is  a minimizer  for  (LP). 

Since  (LP)  has  a feasible  solution  and  [c  — AT/3j]i  ^ 0,  for  each  i e Ib • We  can 
find  a positive  finite  nonzero  aj  at  each  sub-iteration.  This  implies  that 

LZBj(Pi)  < LzBj(Pj+i)  = LzBj+i(j3j+i) 

for  each  sub-iteration  j. 

If  d = 0 at  some  sub-iteration  and  u,-  > u,  or  vt  < /,  for  some  i.  Write  N = 
[B1,  N1], where  B1  is  the  set  of  indices  satisfy  V{  > Ui  or  u,-  < Thus,  by  the 
projection,  we  have 

b — Bzb  — B1zbx  = d + N'vni  . 
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Suppose  d = 0,  then 

b — Bzb  — B1zBi  = Nxvni 

Note  that  d = 0 implies  that  b — Bzb  = Nv.  Hence  Nv  — B1xb i - iV1^  = 0. 
However,  N has  linearly  independent  columns  by  assumption,  it  implies  w,-  = u,  or 
Vi  = li  for  i e /fii,  which  is  contrary  to  the  assumption.  So,  d ^ 0. 

Since  d satisfies  <F(&  - Bzb  — B1zb i)  > 0,  we  have 

dF(b  — Bzb  — Bxzb j)  > 0 
=>■  cF(./Vt>  — B1xbi)  > 0 
=»  <F(A^17  + B1(vbi  - xBi))  > 0 
=P-  drB1(vBi  — xb1)  > 0 

where  u7,  = [ugi,7T]. 

Some  among  those  indices  we  added  in  might  satisfy  (Fb,(vt  - x,)  < 0 which  does 
not  contribute  to  finding  a ’’good”  direction  d.  Hence  we  can  delete  those  indices 
from  IB i.  Note  that  there  is  at  least  one  index  i € IBi  which  satisfies  (Fbt(vi  — Zi)  > 0. 
Without  loss  of  generality,  we  can  just  discuss  the  following  case: 

Suppose  d = 0 at  some  sub-iteration  then  it  follows  that  there  exists  an  index  i 
such  that  Vi  > u;  (similar  for  the  case  when  u,  < /,•).  Let  i be  the  new  index  being 
added  in  the  active  set  from  the  previous  non-active  set  In  and  N = [TV1,  6,].  Again, 
by  the  assumption  that  N has  linearly  independent  columns,  d ^ 0. 

Moreover,  since  (F(b  - Bzb  - biUi)  > 0.  we  have 

<F( b - Bzb  - b^i)  > 0 
=*>  <F(Nv  — b^i ) > 0 
=»  (FN1^  -(-  <F6,(ui  — tt,)  > 0 

=>  d?bi(vi  - Ui ) > 0 
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=>  (F  bi  > 0 

where  vT  = [u,-,7r]. 

Similarly  when  w,-  < /,  we  obtain  dFb{  < 0.  Hence  the  i-th  component  adFbi  of  the 
coefficient  in  the  minimization  part  of  the  Lagrangian  dual  problem  gives  a proper 
sign  with  respect  to  both  cases  ( u,-  > u,  and  u,-  < /,).  This  guarantees  that  a positive 
finite  non-zero  a for  A*+1  will  be  found  which  strictly  increase  L( A).  It  follows  that 

L(\k)  = LzBq( A*)  < Lzg(\j+i)  = L( A*+1)  Vfc. 

hence  Ig  never  repeats  at  each  iteration.  Because  Ig  is  chosen  from  a finite  set,  this 
algorithm  will  terminate  after  finite  many  iterations.  □ 

2.4  Numerical  Experiment 

The  directional  vector  d used  in  each  sub-iteration  of  the  dual  active  set 
algorithm  for  linear  programming  is  determined  by  the  projection  of  the  gradient 
(b  — Bzg ) of  LZB  onto  the  null  space  of  NT . The  formula 

d = (/  - N(NTN)~lNT){b  - Bzb) 

can  only  be  used  under  the  assumption  when  the  columns  of  the  non-active  sub- 
matrix N are  linearly  independent,  otherwise  numerical  errors  may  be  expected. 
Since  we  should  not  check  this  assumption  before  solving  every  problem,  a procedure 
to  separate  the  dependent  columns  from  the  independent  ones  should  be  adopted  in 
the  numerical  implementation.  We  now  discuss  a QR  factorization  approach  in  our 
numerical  implementation. 

Let  Q and  R be  the  QR  factorization  of  a m by  s non-active  sub-matrix  N, 
where  Q = [Qx, . . . , Qm\  is  a m by  m orthonormal  matrix  consists  of  columns  Qi , 
i = 1, . . . , m.  and  R = [Rx, . . . , Rs]  is  a m by  s upper  triangular  matrix  with  columns 
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Ri,  i = l,...,a.  The  diagonal  entries  of  R carry  an  pertinent  information  for  the 
problem.  Let  rjj  denote  the  j- th  entry  on  the  diagonal  of  R.  Since  | ry y | measures  the 
distance  from  column  j of  N to  the  space  spanned  by  the  first  ( j — 1)  columns.  If 
lrii  | = 0,  then  the  first  ( j — 1)  columns  form  a basis  for  the  range  of  N and  the  other 
columns  of  N form  a basis  for  the  null  space  of  NT.  Therefore,  if  Qp  denotes  the 
sub-matrix  consists  of  the  basis  of  the  null  space  of  , then  the  directional  vector 
d can  be  expressed  as 
d = QrQTr(b-BzB). 

We  maintain  a QR  factorization  of  the  non-active  sub-matrix  N at  each  step  and 
update  this  factorization  whenever  a column  is  added  into  N from  the  active  sub- 
matrix B. 

We  implement  the  dual  active  set  algorithm  for  linear  programming  to  solve  45 
test  problems  extracted  from  Netlib.  Netlib  is  a large  electric  software  library  which 
contains  a variety  of  mathematical  software  and  test  problems.  Netlib  files  are  avail- 
able to  public  via  anonymous  ftp. 

The  numerical  results  are  obtained  by  running  program  coded  in  Fortran  on  a SUN 
Spare-station  20  with  optimization  option.  In  this  Fortran  code,  no  attempt  was  made 
to  take  the  advantage  of  the  sparsity  of  the  constraint  matrix  A and  the  coefficient  c. 
Note  that  most  of  the  test  problems  used  in  our  experiment  have  sparsity  less  than 
1%.  Therefore,  comparison  on  the  C.P.U.  time  with  other  numerical  comparison 
results  in,  for  example,  Tomlin  [29]  , is  not  realistic.  However,  from  the  results  in 
table  2.1,  we  are  encouraged  by  the  ratio  of  number  of  total  sub-iteration  and  number 
of  variable,  denoted  by  s/n.  Among  those  test  problems,  34  test  problems  have  this 
ratio  s/n  less  than  2,  11  of  these  34  test  problems  have  even  lower  s/n  ratio  which  is 
less  than  1.  Only  3 test  problems  have  s/n  ratio  larger  than  10.  Analagously,  most 
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of  the  s/m  ratios  (the  ratio  of  the  total  sub-iterations  and  number  of  rows  of  the 
matrix  A)  have  a common  bound. 

We  do  not  take  the  advantage  of  the  sparsity  when  we  code  the  program.  Conse- 
quently this  produces  an  inefficiency  in  the  computation.  With  a restriction  of  both 
hardware  and  software,  we  only  pick  out  those  problems  with  row  number  of  the 
constraint  matrix  A less  than  500.  We  should  expect  to  see  much  better  performance 
leap  for  this  algorithm  provided  a program  coded  with  the  ability  of  processing  the 
sparsity  of  matrix  is  used. 
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Table  2.1.  Performance  data  for  Netlib  test  problems 


Problem 

Name 

m 

n 

Sub- 

iteration^) 

Main 

iterations 

s/m 

s/n 

CPU  Time 
(secs.) 

adlittle 

56 

138 

200 

7 

3.57 

1.45 

.62 

afiro 

27 

51 

59 

4 

2.19 

1.16 

.05 

agg 

488 

615 

832 

16 

1.70 

1.35 

148.04 

bandm 

305 

472 

981 

24 

3.22 

2.08 

89.35 

beaconfd 

173 

295 

396 

6 

2.29 

1.34 

8.45 

blend 

74 

114 

131 

8 

1.77 

1.15 

.63 

bore3 

233 

333 

663 

11 

2.85 

1.99 

28.07 

brandy 

220 

303 

578 

16 

2.63 

1.91 

22.24 

capri 

271 

476 

1887 

20 

6.96 

3.96 

60.31 

degen2 

444 

757 

899 

16 

2.02 

1.19 

221.77 

e226 

223 

472 

860 

25 

3.86 

1.82 

42.20 

etamacro 

400 

734 

2326 

37 

5.82 

3.17 

361.07 

finnis 

497 

1019 

1322 

8 

2.66 

1.30 

255.20 

fitld 

24 

1049 

929 

76 

38.71 

0.89 

7.11 

grow  15 

300 

645 

8008 

53 

26.69 

12.42 

357.46 

grow22 

440 

946 

21765 

100 

49.47 

23.01 

1859.85 

grow7 

140 

301 

917 

10 

6.55 

3.05 

14.02 

israel 

174 

316 

740 

27 

4.25 

2.34 

13.60 

kb2 

43 

68 

80 

7 

1.86 

1.18 

.19 

lotfi 

153 

366 

485 

10 

3.17 

1.33 

9.05 

pilot4 

410 

1153 

14730 

115 

35.93 

12.78 

1126.79 

recipe 

91 

178 

353 

3 

3.88 

1.98 

1.49 

sc  105 

105 

163 

118 

7 

1.12 

0.72 

1.05 

sc205 

205 

317 

211 

7 

1.03 

0.67 

7.09 

sc50a 

50 

78 

53 

2 

1.06 

0.68 

.12 

sc50b 

50 

78 

55 

2 

1.10 

0.71 

.12 

scagr25 

471 

671 

1144 

47 

2.43 

1.70 

245.84 

scagr7 

129 

185 

263 

44 

2.04 

1.42 

4.27 

scfxml 

330 

600 

1800 

32 

5.45 

3.00 

163.04 

scorpion 

388 

466 

452 

4 

1.16 

0.97 

52.34 

scrs8 

490 

1275 

2588 

39 

5.28 

2.03 

578.26 

scsdl 

77 

760 

152 

2 

1.97 

0.20 

1.25 

scsd6 

147 

1350 

646 

9 

4.39 

0.48 

15.95 

scsd8 

397 

2750 

3736 

21 

9.16 

1.32 

453.22 

sctapl 

300 

660 

834 

10 

2.78 

1.26 

45.79 

sharelb 

117 

253 

353 

16 

3.02 

1.40 

7.98 

share2b 

96 

162 

318 

18 

3.31 

1.96 

2.14 

ship041 

402 

2166 

509 

6 

1.27 

0.23 

267.00 

ship04s 

402 

1506 

616 

8 

1.53 

0.41 

191.07 

stair 

356 

534 

697 

13 

1.96 

1.31 

72.70 

standata 

359 

1258 

418 

2 

1.16 

0.33 

85.74 

standmps 

467 

1258 

1340 

6 

2.87 

1.07 

261.90 

stocforl 

117 

165 

182 

11 

1.56 

1.10 

2.07 

tuff 

333 

627 

1387 

17 

4.17 

2.21 

92.20 

vtp.base 

198 

329 

519 

11 

2.62 

1.58 

14.75 

CHAPTER  3 

CONJUGATE  GRADIENT  ACTIVE  SET  ALGORITHM 
3.1  Introduction 

Quadratic  programming  is  a fundamental  problem  in  optimization  and  in 
general  can  be  used  as  basic  subroutine  for  the  nonlinear  optimization  algorithm.  For 
example,  see  Conn,  Gould,  and  Toine  [4].  In  this  chapter,  we  present  an  algorithm 
for  solving  a specific  quadratic  programming  problem  with  box  constraints  which  can 
be  expressed  as  (QP): 

min  f{x ) = \)xtQx  + qTx  (3.1) 

Li 

s.i.  I < x < u,  x 6 Rn ■ 

where  Q is  a symmetric  positive  definite  matrix  and  /,  u , and  q € Rn.  This  algorithm, 
namely,  the  conjugate  gradient  active  set  algorithm,  combines  the  primal  active  set 
strategy  with  the  finite  convergence  property  of  the  conjugate  gradient  method  for 
quadratic  programming. 

We  introduce  the  conjugate  gradient  active  set  algorithm  in  section  3.2.  A modi- 
fied version  of  this  algorithm  in  the  procedure  to  update  the  active  indices  is  presented 
too.  We  also  point  out  the  relation  between  this  algorithm  and  the  primal  active  set 
algorithm  presented  in  section  1.2.  In  section  3.3,  we  show  that  the  conjugate  gradi- 
ent active  set  algorithm  terminates  at  an  optimal  solution  of  the  quadratic  problem 
(QP)  in  finitely  many  iterations.  In  the  last  section,  we  compare  the  relative  speed 
of  the  conjugate  gradient  active  set  algorithm  and  the  following  algorithms  on  a test 
problem  of  the  form  (QP)  from  engineering: 

• Gradient  projection  algorithm  of  Goldstein  [11]  and  Levitin  and  Polyak  [22]. 

• Gradient  Projection  algorithm  of  Rosen  [26],  [27]. 
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• Interior  point  algorithm  of  Han,  Pardalos,  and  Ye  [19]. 

A collection  of  8 different  problems  are  obtained  by  varying  both  the  box  constraints 
and  the  Hessian  matrix  of  the  original  test  problem.  A brief  introduction  of  each 
algorithm  tested  in  this  section  is  followed  by  a report  of  the  numerical  experiment 
results  and  figures  for  the  optimal  solution  of  each  test  problem. 

3.2  The  Algorithm 

In  this  section,  we  present  an  algorithm,  namely,  the  conjugate  gradient  active 
set  algorithm,  to  minimize  a quadratic  programming  problem  with  box  constraints. 
We  consider 


min  \xTQx  + qTx  (3.2) 

s.t.  I < x < u 

where  Q is  a n by  n positive  definite  matrix  and  /,  u,  and  q are  constant  vector  in 
Rn . The  conjugate  gradient  active  set  algorithm  is  an  implementation  of  the  primal 
active  set  algorithm  presented  in  section  1.2.  The  sub-iteration  procedures  (a)  and 
(b)  in  the  primal  active  set  algorithm  are  assured  by  the  conjugate  gradient  iterations 
involved  for  solving  an  unconstrained  quadratic  problem  at  each  sub-iteration. 

We  introduce  notations  and  an  auxiliary  problem  at  first. 

Notation: 

1.  The  index  set  Ib  and  In  and  the  vector  xb  are  defined  the  same  as  in  chapter 

2. 


2.  Given  a vector  x,  define  Pb  be  the  projection  that 


o ; i € B 


[Pbx]x  = < 


Xi  ; otherwise 
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3.  Given  an  index  set  /#,  define  an  auxiliary  problem  as  follow: 

min  ±xTQx  + qTx  (3.3) 

s.t.  x b = zb 

where  x,-,  i 6 N , is  unconstrained. 

4.  Denote  <7,-  = g(z{)  as  the  gradient  of  the  cost  function  evaluated  at  Z{. 

The  conjugate  gradient  active  set  algorithm  starts  at  a feasible  point.  An  active 
index  set  is  determined  by  the  collection  of  all  the  entries  of  the  initial  point  which 
are  at  bound.  An  unconstrained  quadratic  subproblem  (SQP)  is  solved  according  to 
the  active  index  set  by  the  conjugate  gradient  method.  We  use  the  Polak-Riebiere 
form  for  the  conjugate  gradient  method  in  this  algorithm. 

The  Polak-Ribiere  form  of  the  conjugate- gradient  algorithm  applied  to  problem 
(SQP)  can  be  implemented  in  the  following  way: 


Polak-Ribiere  form 

Set  d0  — -PbS'o,  and  at  iteration  i,  we  have 


Zi+l 


Cti 


di+ 1 


Pi 


Zi  -(-  a,  d,  , 

di9i 

W 

— Pb5»+ 1 + M, 
gf+iPB(g>+i  - g,) 
gJPBgi 


(3.4) 

(3.5) 

(3.6) 

(3.7) 


If  any  solution  of  the  quadratic  subproblem  is  going  out  of  the  box  constraint, 
we  extend  the  active  index  set  by  choosing  the  step  size  along  the  directional  vector 
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until  it  first  reaches  the  boundary  of  the  box  constraints.  This  is  the  procedure  to 
update  (add  in)  the  active  index  set.  Also,  if  at  some  iteration,  the  conjugate  gradient 
iteration  attains  a solution  x * of  a current  quadratic  subproblem  (SQP)  which  is  also 
feasible  in  the  box  constraints  / < x*  < u,  a procedure  is  used  to  update  (drop  out) 
the  active  index  set. 

We  present  the  conjugate  gradient  active  set  algorithm  as  follow: 


Conjugate  Gradient  Active  Set  Algorithm 

Given  initial  guess  x0,  set  dQ  = —Psg{x 0). 

Main  iteration: 

Given  Xky  let  j = 0,  yj  = z0  = Xk,  and  define 

Ib0  = {*  : (^o).'  = U and  g(z0)  > 0,  or  ( z0)i  = u,-  and  g(z0)  < o} 

Sub-iteration: 

Apply  the  Polak-Ribiere  form  of  the  Conjugate  Gradient  method  to  (SQP) 

(i) .  If  at  some  iteration,  z,-+1  = z,-  + a,d,  is  infeasible  in  problem  (QP),then  let  a be 

the  largest  a with  the  property  that  z,  + ad,  is  feasible,  set  z = z,-  + ad,-,  let 
Ibj+ i be  the  set  defined  by 

Ib}+1  = IB,  U {s  : z.  = l,  or  z,  = u,}, 

set  yj+i  = z,  increment  j , and  proceed  the  sub-iteration. 

(ii) .  If  at  some  iteration,  z,+1  is  a solution  to  problem  (SQP)  or  number  of  the 

Conjugate  Gradient  iterations  exceeds  number  of  variables  in  problem  (SQP), 
then  set  xjt+i  = z,+i,  and  proceed  to  the  next  iteration. 
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We  now  modify  the  conjugate  gradient  active  set  algorithm  so  that  the  produre 
of  adding  in  indices  to  the  active  index  set  is  more  efficient.  When  a solution  of  the 
quadratic  subproblem  (SQP)  is  going  out  of  the  box  constraints,  a line  search  routine 
is  employed  to  determine  a suitable  step  size  along  the  direction,  then,  a projection 
of  this  new  point  into  the  box  constraints  enable  the  algorithm  to  add  in  indices  to 
the  active  index  set. 

Hager  [14]  discussed  a combined  Armijo/ Goldstein  line  search  routine  which  is 
adopted  in  the  line  search  for  the  gradient  projection  method  of  Goldstein  et  al. 
in  this  numerical  experiment  and  stated  in  section  3.4.1.  We  state  the  modified 
conjugate  gradient  active  set  algorithm  as  follow: 


Modified  Conjugate  Gradient  Active  Set  Algorithm 
Modify  procedure  (i)  in  the  original  algorithm  to  (i’). 

(P).  If  at  some  iteration,  z;  + a,d,  is  infeasible  in  problem  (QP),then  let  a be  the 
largest  a with  the  property  that  z,-  + ad,  is  feasible,  set  z = z,-  + ad,. 

Denote  z,+i  = P [z,-  + s(z  — (z,-  + a,d,)] , 

where  P is  a projection.  Use  the  Armijo/Goldstein  line  search  routine  to  de- 
termine a step  size  s. 


The  projection,  denoted  as  P used  here  is  described  as  follow: 

Given  a vector  x , the  projection  of  x;  into  the  box  constraints  is  Px  where 


[Px],-  = /,  if  x,  < /,• 

[Px]t-  = x,-  if  /,  < Xi  <Ui 
[Px]f-  = Ui  if  x,  > U(. 
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This  procedure  may  add  in  more  than  one  index  to  the  active  index  set  in  one 
iteration.  Numerical  experiment  results  also  show  that  it  speed  up  the  convergence 
rate  of  the  conjugate  gradient  active  set  algorithm. 

3.3  Convergence 

We  can  show  that  the  conjugate  gradient  active  set  algorithm  attains  an 
optimal  solution  to  the  quadratic  programming  problem  (QP)  in  a finite  number  of 
iterations.  In  the  proof  of  the  following  theorem,  we  denote  Igk  as  the  last  active 
index  set  of  the  fc-th  main  iteration. 

Theorem  3.3.1  The  conjugate  gradient  active  set  algorithm  for  the  quadratic  problem 
3.3  will  terminate  after  finitely  many  iterations  and  the  generated  sequence  {a;*}  will 
end  at  a point  x*  such  that  V f(x*)T(x  — x*)  > 0 for  all  x satisfies  l < x < u. 

Proof  of  Theorem  3.3.1  We  first  show  that  the  sub-iteration  process  under  each  main 
iteration  will  be  terminated  after  finitely  many  iterations,  and  the  objective  function 
is  strictly  decreased  after  each  main  iteration. 

At  the  k- th  iteration,  in  the  sub-iteration  process  (i),  if  Z;+i  = 2, -fad,  is  infeasible 
in  problem  (3.1)  and  a is  the  largest  step  size  such  that  z = z,-  + ad{  is  feasible,  then 
0 < a < a.  Since  /(x)  is  strictly  convex  and  /(z,+ 1)  < /(z,),  we  have 

m = /((i--)2,+^2l+1) 

a a 

< (1  - i)f{Zi)  + f/(zi+1) 

a a 

< (1  - §)/(*)  + £ /(*•■) 

a a 

= /fa)- 

Therefore  /(z)  < /(z,).  Moreover,  the  set  Ib3  is  strictly  contained  in  7bj+1.  Since 
Jb>+1  is  a subset  of  {1,2,  ...,n},  this  sub-iteration  process  will  be  terminated  after 
finitely  many  iterations. 
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In  the  sub-iteration  process  (ii),  if  z,  = z terminates  this  sub-iteration  process 
then  [v/(^)]«'  = 0 for  i 6 /#*.  By  Taylor’s  Theorem  and  the  fact  that  Q is  a positive 
definite  matrix,  for  any  i £ if,  i / z,  and  xgk  = zgk,  we  can  write  x = z + d where 
dgk  = 0 and  dfjk  ^ 0.  Then  there  exists  a positive  A E (0, 1)  such  that 

/(*)  ~ f(z)  = \/f(z)Td  + A dTQd 

— f(z)]J}kdftk  + [syf(z)]gkdgk  + A dT  Qd 
= A dTQd 
> 0. 

hence,  z,  is  a strict  minimizer  for  (ii).  Moreover,  the  Conjugate  Gradient  method 
ensures  that  this  sub-iteration  will  stop  after  finitely  many  iterations. 

Since  the  function  value  of  f(x)  is  strictly  reduced  after  each  main  iteration,  B ^ 
never  repeats.  Therefore,  by  the  finitely  many  choices  of  B W from  the  set  {1,2,  .. , n}, 
the  conjugate  gradient  active  set  algorithm  will  reach  a minimizer  x*  after  finitely 
many  iterations. 

If  this  algorithm  has  not  terminated  at  the  k- th  iteration,  then  /B*+i  must  be 
strictly  contained  in  Igk,  otherwise  /B*+i  = Igk.  If  /B*+i  = Igk,  then  for  all  x 
satisfies  / < x < u 

V/(xfc+1)T(x  - x*+1) 

= [V/(3:<!+1)]b(*+i)(2;  — Xk+1)B(k+i)  + [V/(x*+1)]]vB(*+i)(iC  — Xk+1)Ng(,k+i) 

= [Vf(xk+1)]TgW(x  - xk+1)BW  + [vf(xk+1)]ZgW(x-xk+')NgW 
= [V/(x*+1)]fi(*+i)(x  - x*+1)j5(*+i) 

> 0. 

This  implies  that  \7f(xk+1)T(x  - xk+1)  > 0 for  all  x satisfies  l < x < u.  □ 
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3.4  Numerical  Experiment 

An  interesting  application  of  quadratic  programming  with  box  constraints  is 
the  obstacle  problem  which  is  to  find  the  equilibrium  position  of  an  elastic  membrane 
subject  to  a vertical  force  and  an  obstacle  bounded  above  and  below.  This  problem 
can  be  formulated  as 

minf(x)  = i xtQx  + qTx  s.t.  I < x < u,  (3.8) 

where  x,  u,  l,  q £ Rn , and  Q £ RnXn  is  a structured  block  tridiagonal  matrix  which 
is  symmetric  positive  definite.  All  the  entries  off  the  main-diagonal,  super-diagonal, 
and  sub-diagonal  block  of  this  matrix  are  zero. 


where 


Q = 


T = 


T -I 
-I  T -I 
-I 


T -I 


4 -1 
-1  4 -1 

-1 


-/  T -I 
-I  T 


4 -1 


-1 


4 -1 
-1  4 


The  I matrix  in  Q is  a m x m identity  matrix,  T € ir*xm,  and  m2  = n.  The 


vector  q in  f(x)  is  given  by  q = -h2  where  h = 

In  the  computational  experiment,  we  also  test  three  more  problems  each  with 
the  Q matrix  perturbed  by  multiplying  a diagonal  matrix  D to  both  sides  of  Q. 
Accordingly  we  have  the  following  three  cases  for  the  matrix  D: 


1.  The  entry  at  (1,1)  is  100,  the  rest  entries  of  the  diagonal  are  all  1. 


43 


2.  The  matrix  D is  formulated  into  a block  diagonal  matrix,  each  block  matrix  is 
a m x m diagonal  matrix  formulated  by 

1 


(*»*')  = 1 - 


2 x (m  + 1 - i) 


7T  i = 1, . • . ,m. 


3.  The  diagonal  entries  of  D are  randomly  generated  from  (0, 1). 


Two  different  bound  constraints  are  considered  for  each  of  these  four  quadratic  prob- 
lems. 

I: 


Ui  = (sm(9.2a,)  x s*n(9.37,-))2  + 0-02 
/,  = (sm(9.2a,)  x sin(9.37,))3, 


Ui  = 2000 

/,•  = sm(3.2a,)  x sm(3.37,), 


where 


ii  — 1) 

a*  = (* x m)  x h, 


m 


i x h 

7 * = , 


m 


for  t = 1, . . . ,n. 

Therefore,  we  have  total  of  four  test  problems  and  two  sets  of  constraints  for 
each  test  problem.  For  convenience  of  notation,  we  denote  these  problems  by  prob- 
lem (constraint  set  number)-(problem  number).  For  example,  problem  II-3  denotes 
problem  3 with  the  second  constraint  set. 

The  algorithms  included  in  this  experiment  are 
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• Conjugate  gradient  active  set  method. 

• Conjugate  gradient  active  set  method  + Linear  Search. 

• Gradient  projection  method  of  Rosen. 

• Gradient  projection  method  Goldstein  et  al. 

• Interior  point  method  of  Han  et  al. 

In  order  to  make  the  starting  point  feasible  in  the  box  constraints,  we  set  the 
starting  point  x°  = |(/  + u)  in  every  program.  Therefore  the  starting  point  is  feasible 
at  the  first  iteration.  We  also  use  the  same  stopping  criterion  in  the  program  for  each 
algorithm.  The  stopping  criterion  is  considered  in  the  following  way: 

At  a local  minimizer  x*  of  the  problem 

minf(x)  s.t.  I < x <u  (3.9) 

we  have 

1.  If  /,  < x*  < Ui,  then  ^ = 0 

OXi 

2.  If  x*  = then  > 0 

OXi 

3.  If  x*  = U{,  then  --}X  ^ < 0 

OXi 

Our  error  expression  measures  how  much  these  conditions  are  violated  at  x where 
x is  feasible.  Therefore  the  tolerance  function  E(x)  is  given  by 

1 

2 

(3.10) 


E(x)  = 


E ( 


dXi  J 


where 
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N = 

{*  • It  < 

X;  < it,-} 

A , = 

{* : = 

Au  = 

{*  : Xi  = 

u,} 

Aj-  = 

{ it  A 

T<°> 

K = 

{*  € Au  : 

T>«>- 

All  test  programs  will  be  stopped  when  the  value  of  the  tolerance  function  E(x) 
reaches  10-4. 

Among  the  algorithms  we  test  in  this  numerical  experiment,  the  interior  point 
method  for  solving  the  quadratic  programming  problem  with  box  constraints  can  be 
found  in  [19]  by  Han,  Pardalos,  and  Ye.  They  also  implemented  this  method  on 
solving  some  engineering  problems  by  utilizing  the  vector  processing  ability  on  an 
IBM3090  — 600S.  The  code  we  obtained  from  H-D  Chen  is  a rewritten  version  of 
their  code  to  work  on  our  SUN  workstation. 

The  gradient  projection  method  of  Rosen  [26]  [27]  , projects  the  negative  gradient 
to  reduce  the  cost  function  and  still  keep  the  new  iteration  point  inside  the  feasible 
region.  Each  iteration  is  obtained  by 

xk+1  = xk  — dkPgk 

where  the  projection  matrix  P is  determined  by  the  active  constraints  at  each  iter- 
ation. The  step  size  a*  is  chosen  to  minimize  the  cost  function  along  the  direction 
—Pgk  and  also  maintain  the  feasibility  at  the  same  time. 

The  idea  of  this  gradient  projection  method  has  been  widely  used  for  solving 
nonlinear  programming  problems.  Despite  being  one  of  the  most  practical  primal 
methods,  the  convergence  of  this  method  had  been  remained  open  for  more  than  25 
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years  until  a complete  proof  of  Du  and  Zhang  [7]  appeared  in  1989.  They  established 
the  convergence  result  by  showing  that  a non-KKT  point  cannot  become  a limit  point 
of  the  sequence  generated  by  the  gradient  projection  method. 

In  our  test  problem  the  projection  matrix  P,  due  to  the  simple  structure  of  the 
box  constraints,  is  a diagonal  matrix.  The  value  of  each  entry  p„  on  the  diagonal  of 
the  matrix  P is 


0 


Pit  = 


Xi  — Ui 


Or  Xi  = 1 


1 ; otherwise. 

Exact  line  search  is  used  for  adding  active  indices  in  this  algorithm.  We  drop  out 
indices  from  the  active  index  set  which  satisfy  one  of  the  following  rules: 


1.  x(i)  = u(i),  g{i)  > 0,  and\g(i)\  > 10  • max\g(j)\.  (3.11) 

2.  x(z)  = l(i),  g(i)  < 0,  and  |sr(i)|  > 10  • max\g(j)\.  (3.12) 

In  the  next  sub-section  we  briefly  describe  the  gradient  projection  method  of 
Goldstein  et  al.  and  the  line  search  routine  we  used  in  this  numerical  experiment. 

3.4.1  Gradient  Projection  Method  of  Goldstein  et  al.  and  Line  Search 

The  gradient  projection  method  of  Goldstein  [11]  1964  and  Levitin  and  Polak 
[22]  1965  consists  of  the  following  iterations  : 

xk+1  = P(xk  - ak9k)  (3.13) 

where  ak  is  step  size  and  the  projection  P in  (3.13)  is  defined  to  project  a point  onto 
the  feasible  region  by 


P(Xi)  = 


Xi,  if  li  < Xi  < Ui 
li,  if  Xi  < /,• 

Ui,  if  Xi  > Ui 
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Numerous  articles  concerning  the  convergence  properties  of  this  gradient  projec- 
tion method  have  been  published.  Bertsekas  established  a finite  terminated  gradient 
projection  method  based  on  an  Armijo  type  line  search  routine.  An  important  re- 
sult of  Bertsekas’s  work  is  that,  if  a sequence  of  points  generated  by  this  method 
converges  to  a local  minimizer  which  satisfies  the  strict  complementarity  and  second 
order  sufficiency  conditions,  then  the  active  constraints  at  the  local  minimizer  can  be 
identified  in  a finite  number  of  iterations.  Calamai  and  More  [3]  1987,  had  the  same 
result  on  a linearly  constrained  problem  provided  the  gradient  projection  method 
converges  to  a non- degenerate  KKT  point. 

We  discuss  the  line  search  routine  as  follows: 

Let  f(x)  be  a convex  C2  function.  Consider  the  problem 

minf(x)  s.t.  l<x<u  (3.14) 

by  using  the  projected  gradient  method  above. 


Line  Search  (I) 

The  step  size  ak  in  (3.13)  is  chosen  by  setting  ak  = /3ma  where  a £ (0,  oo),  /?  £ (0, 1) 
and  m £ Z.  At  the  (k+l)th  iteration,  we  have  xk+i  determined  by  (3.13),  where  the 
a in  ak  is  initial  set  to  be  equal  to  ak-\  and  m = 0 for  the  initial  value  of  ak.  then 
adjust  m by  following  the  rule. 

Fix  c £ (0, 1), 

Consider 

/(**)  - /(* fc+i)  > e • V/(**)T  • (xk  ~ Xfc+i)  (3.15) 

If  the  inequality  (3.15)  is  true,  then  we  have  either  (a)  or  (6): 
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(a)  decrease  m from  —1  until  the  inequality  first  fails  at  m = — p.  use  m = 

— (p  — i);  stop. 

(b)  decrease  m until  each  component  of  xk+i  reaches  the  bound. i.e.  xj^  = «(*) 

or  = /(')  for  each  i. 

Otherwise,  increase  m from  1 until  it  satisfies  the  inequality. 


Proposition  3-4-1  The  rule  for  choosing  ak  is  well  defined  and  a*  can  be  obtained  in 
finitely  many  computations  of  the  inequality  (3.15). 

Proof  of  Proposition  3.4-1  If  the  inequality  is  true  and  m = — (p  — 1)  for  some  p. 
then  it  stops  in  finitely  many  of  function  evaluations,  otherwise  it  will  stop  when 


each  component  of  xk+i  reaches  the  bound.  In  this  case  m is  finite. 

On  the  other  hand,  If  the  inequality  fails, that  is,  we  have 

/(**+ 1)  > f(xk ) + e • V/(z*)r  • (*M-i  - xk)  (3.16) 

then  since  /(x)  is  a C 2 function,  we  have 

/(**+ 1)  = /(**)  + V/(**)T  • (**+i  ~ xk)  + 0(\xk+1  - x*|2)  (3.17) 

Note  that  (3.16)  and  (3.17)  imply  that 

(1  - e)  * Vf(*k)T  • (**+1  - *k)  + 0{ \xk+1  - xk\2)  > 0 (3.18) 

Also  by  Note  1 (next  page)  and  0 < t < 1,  we  have 

(1  - e)  • V/(**)r  • (**+i  - xk)  < 0 (3.19) 


For  sufficiently  large  m such  that  0(|x*+1  - xk\2)  = 0,  inequality  (3.18)  and  (3.19) 
give  a contradiction.  Therefore,  ak  must  be  obtained  in  finite  many  computations  of 
the  inequality.  □ 
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Theorem  3-4.1  If  {x,}  is  a sequence  of  points  generated  by  (3.13)  and  each  ak  is 
chosen  by  the  line  search  (I),  then  every  limit  point  of  {xk}  is  a stationary  point. 

Proof  of  Theorem  3.1.1  If  {x,}  converges  to  a stationary  point  we  are  done.  Hence 
assume  otherwise.  Suppose  that  {x,}  is  a sequence  of  points  generated  by  the  above 
rule  and  {xk}k€i  is  a subsequence  converging  to  a point  x*.  Since  {/(x*)}  is  a 
monotonic  decreasing  sequence.  f(xk)  -*  /(x*)  and  \f(xk)  - /(x*+1)|  — >•  0.  For  any 
y,  l < V < u-  we  have 

V/(**)r  • (x*  - y) 

= v/(x*;)r  • (xk  - Xfc+i)  + y/(x k)T  • (**+1  - y) 

< S/f(xk)T  ■ ( Xk  - xk+1)  + — • (xjfe  - xk+1)T  ■ (®*+1  - y) 

ak 

< — - + v/(x*)l 

ak 

Let  k -4  oo,  we  have  xk  ->  x * and  \xk  - ®fc+1|  ->  0.  hence  y f(**)T  • (x*  - y)  < 0 for 
any  y in  the  feasible  domain.  □ 

Note: 

1.  For  any  xk  , y in  the  feasible  domain,  define  xk+\  as  (3.13).  Then 

ak  • V/(xfc)T  • [v  ~ z*+i]  > [xk  ~ Zfc+i]7,  • [y  - xjt+i] 

2. 

f(xk)-f(xk+ 1)  > C-  v/(xjt)r  • (xfc  -Xk+l) 

> e.  l£*z£*±il! 

— ak 

We  employ  a combined  Armijo/Goldstein  line  search  routine  proposed  by  Hager 
[14]  in  the  gradient  projection  method  of  Goldstein  et  al.  in  the  numerical  ex- 
periment. We  state  Armijo’s  rule  and  Goldstein’s  rule  followed  by  the  combined 
Armijo/Goldstein  line  search  routine.  Denote  x(a)  = P(x  — a y f(x)). 
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Armijo’s  rule:  Given  7 > 1,  0 < e < 1,  find  a step  size  a satisfies  both 


f(xk)  ~ f(xk(a))  > e V f{xk)(xk  - xk(a)), 


and 


f(xk)  - /(**( 7«))  > e V f(xk)(xk  ~ xk(ja)). 


Goldstein’s  rule:  Given  0 < p < 0.5,  find  a step  size  a satisfies 


(!  ~ P)  V f(xk)(xk  - xk(a))  > f(xk)  - f(xk(a))  > p y f(xk)(xk  - xk(a)). 


Line  Search  (II) 

Given  a starting  step  size  a equal  to  the  previous  step  size  and  a positive  integer 
M,  let  ak  = (3ma , where  m is  determined  by 

(a) .  If  f(xk)  — f(xk(a))  < e V f{xk)(xk  — xk(a)),  then  choose  m as  the  first  integer 

< 0 such  that  the  Armijo’s  rule  hold. 

(b) .  If  f(xk)  — f(xk(/3ma))  > e V f(xk)(xk  — xk(/3ma)),  then  starts  at  m = 0,  find 

the  first  nonnegative  integer  m such  that 

(1  - P)  V f(xk)(xk  ~ xk(r<*))  > f(xk)  - f{xk[pra)). 

If  m > M,  then  set  ak  = (3Ma. 

lim  < M and  f(xk)  - f(xk(/3ma))  > p y f(xk)(xk  - xk(/3ma)),  then  ak  = (3ma 
satisfies  the  Goldstein’s  rule. 

Otherwise,  ak  = /?m-1a  satisfies  the  Armijo’s  rule. 
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3.4.2  Numerical  Conclusions 

The  results  of  this  numerical  comparison  are  obtained  from  a SUN  Sparc  20 
workstation.  The  C.P.U.  run  time  for  each  algorithm  on  all  the  test  problems  appear 
in  table  3.1  through  table  3.8.  The  figure  3.1  and  figure  3.2  provide  three-dimensional 
graphs  for  the  optimal  solution  for  each  of  the  test  problem. 

From  the  C.P.U.  run  time,  we  observe  that  the  conjugate  gradient  active  set 
algorithm  with  Armijo/ Goldstein  line  search  routine  and  Goldstein  et  al.  projection 
performs  very  well  on  most  test  problems.  The  interior  point  method  is  relatively 
slow  on  most  test  problems.  The  C.P.U.  time  used  on  each  problems  are  very  close. 
We,  however,  think  that  other  methods  for  solving  the  linear  systems  may  reduce 
the  C.P.U.  time  as  compared  to  the  Gauss-Seidel  method  employed  in  our  code.  The 
gradient  projection  method  of  Goldstein  et  al.  with  Armijo/Goldstein  line  search 
routine  provides  a mixture  result.  It  performs  extremely  fast  on  some  test  problems, 
however,  it  can  be  affected  by  the  ill-condition  of  the  Hessian  matrix  of  the  respective 
problem.  The  performance  of  the  gradient  projection  method  of  Rosen  is  also  affected 
by  the  ill-condition  of  the  test  problems. 


52 


Tab 

e 3.1.  C.P.U.  time  for  the  test  problem 

1-1 

m 

Rosen 
Gradient  Proj. 

Interior 

Point 

Goldstein  et  al. 
Gradient  Proj. 

Conjugate 

Gradient 

Conjugate  Grad. 
+ Line  Search 

10 

0.03 

0.06 

0.00 

0.02 

0.01 

20 

0.23 

0.82 

0.03 

0.17 

0.05 

30 

1.26 

4.20 

0.15 

0.79 

0.18 

40 

3.61 

12.45 

0.39 

2.17 

0.50 

50 

7.75 

28.83 

0.88 

4.89 

1.00 

60 

16.83 

62.43 

1.51 

10.10 

1.65 

70 

28.37 

115.60 

2.47 

17.67 

2.61 

80 

47.23 

191.72 

3.82 

31.25 

4.07 

90 

78.49 

301.86 

5.33 

50.41 

5.30 

100 

117.03 

485.82 

7.18 

73.93 

9.22 

Tab 

e 3.2.  C.P.U.  time  for  the  test  problem 

1-2 

m 

Rosen 
Gradient  Proj. 

Interior 

Point 

Goldstein  et  al. 
Gradient  Proj. 

Conjugate 

Gradient 

Conjugate  Grad. 
+ Line  Search 

10 

.02 

0.07 

0.01 

0.02 

0.02 

20 

.32 

1.00 

0.04 

0.22 

0.07 

30 

1.66 

5.23 

0.19 

1.04 

0.22 

40 

1363.58 

15.44 

2706.37 

2.95 

0.64 

50 

4777.09 

35.82 

8866.32 

6.43 

1.72 

60 

79.87 

16175.16 

13.11 

9.85 

70 

144.23 

22.38 

5.53 

80 

240.05 

38.57 

5.15 

90 

406.54 

62.08 

7.61 

100 

588.97 

91.29 

9.85 
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Tab 

e 3.3.  C.P.U.  time  for  the  test  problem 

1-3 

m 

Rosen 
Gradient  Proj. 

Interior 

Point 

Goldstein  et  al. 
Gradient  Proj. 

Conjugate 

Gradient 

Conjugate  Grad. 
+ Line  Search 

10 

0.03 

0.08 

0.01 

0.03 

0.01 

20 

0.31 

0.90 

0.04 

0.23 

0.07 

30 

1.48 

4.60 

0.19 

1.05 

0.22 

40 

4.21 

16.40 

0.48 

2.97 

0.63 

50 

10.01 

37.94 

1.05 

6.60 

1.44 

60 

20.65 

75.75 

1.87 

13.04 

2.19 

70 

36.06 

137.18 

3.05 

23.52 

3.92 

80 

57.72 

228.04 

4.63 

39.56 

5.73 

90 

89.32 

393.98 

6.38 

62.46 

7.64 

100 

133.39 

586.66 

8.60 

93.10 

11.76 

Tab 

e 3.4.  C.P.U.  time  for  the  test  problem 

1-4 

al 

m 

Rosen 
Gradient  Proj. 

Interior 

Point 

Goldstein  et  al. 
Gradient  Proj. 

Conjugate 

Gradient 

Conjugate  Grad, 
-f  Line  Search 

10 

0.03 

0.08 

0.03 

0.03 

0.03 

20 

.75 

0.61 

0.64 

0.32 

0.26 

30 

2.53 

3.78 

1.51 

1.15 

0.74 

40 

16.16 

10.53 

2.94 

4.38 

2.64 

50 

11.59 

19.63 

2.15 

8.02 

5.48 

60 

49.33 

29.32 

1.73 

19.25 

10.35 

70 

36.06 

43.23 

2.27 

30.00 

17.73 

80 

68.45 

96.86 

2.86 

59.81 

24.62 

90 

92.67 

119.07 

3.56 

66.72 

43.00 

100 

145.03 

152.44 

4.37 

144.05 

53.79 
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Table  3.5.  C.P.U.  time  for  the  test  problem 


m 

Rosen 
Gradient  Proj. 

Interior 

Point 

Goldstein  et  al. 
Gradient  Proj. 

Conjugate 

Gradient 

Conjugate  Grad. 
+ Line  Search 

10 

0.06 

0.36 

0.03 

0.02 

0.02 

20 

1.09 

4.08 

0.50 

0.42 

0.23 

30 

5.98 

20.42 

2.65 

2.05 

1.40 

40 

19.15 

61.80 

8.56 

7.18 

3.66 

50 

47.56 

150.65 

21.08 

16.46 

8.18 

60 

100.70 

307.25 

44.00 

32.93 

18.88 

70 

188.67 

555.87 

82.23 

59.28 

25.60 

80 

321.31 

957.33 

140.26 

100.05 

43.01 

90 

512.43 

1494.04 

226.25 

165.60 

67.48 

100 

777.50 

2254.83 

347.67 

250.39 

105.52 

1-1 


Table  3.6.  C.P.U.  time  for  the  test  problem  II-2 


m 

Rosen 
Gradient  Proj. 

Interior 

Point 

Goldstein  et  . 
Gradient  Proj. 

Conjugate 

Gradient 

Conjugate  Grad. 
+ Line  Search 

10 

11.16 

0.62 

87.63 

0.03 

0.03 

20 

98.79 

5.33 

1671.68 

0.54 

0.28 

30 

343.18 

25.05 

8980.68 

2.82 

1.67 

40 

825.88 

74.55 

7.94 

4.23 

50 

1705.00 

177.13 

19.41 

9.98 

60 

365.59 

38.53 

20.96 

70 

662.38 

69.37 

29.52 

80 

1108.81 

131.15 

58.51 

90 

1765.15 

204.34 

95.75 

100 

2729.66 

294.25 

111.13 
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Table  3.7.  C.P.U.  time  for  the  test  problem 


m 

Rosen 
Gradient  Proj. 

Interior 

Point 

Goldstein  et  al. 
Gradient  Proj. 

Conjugate 

Gradient 

Conjugate  Grad. 
+ Line  Search 

10 

0.09 

0.45 

0.04 

0.04 

0.04 

20 

1.43 

4.94 

0.68 

0.50 

0.38 

30 

7.24 

25.69 

3.51 

2.65 

1.65 

40 

22.91 

75.76 

11.00 

8.02 

4.65 

50 

55.33 

184.99 

26.86 

19.15 

10.69 

60 

114.20 

375.91 

55.46 

35.36 

20.64 

70 

210.33 

694.18 

105.56 

71.35 

36.64 

80 

371.31 

1150.26 

184.18 

125.64 

54.57 

90 

612.66 

1856.84 

296.23 

206.41 

81.60 

100 

938.27 

2770.96 

452.03 

318.73 

134.43 

1-3 


m 

Rosen 
Gradient  Proj. 

Interior 

Point 

Goldstein  et  al. 
Gradient  Proj. 

Conjugate 

Gradient 

Conjugate  Grad. 
+ Line  Search 

10 

3.32 

0.44 

2.81 

0.15 

0.17 

20 

744.43 

6.36 

388.77 

8.22 

9.63 

30 

1839.16 

38.11 

963.34 

40.02 

39.97 

40 

9131.43 

120.02 

149.38 

116.25 

50 

368.29 

380.79 

379.74 

60 

688.76 

503.70 

461.30 

70 

1665.89 

1554.12 

1319.79 

80 

2467.68 

1816.15 

1513.33 

90 

4552.76 

2868.41 

2602.75 

100 

5179.64 

4168.94 

3418.86 

1-4 
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Figure  3.1.  3-D  graphs  for  the  optimal  solution  to  the  test  problems  with  the  first 
set  of  box  constrains 
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0 0 0 0 
Solution  for  problem  11-3  Solution  for  problem  11-4 

Figure  3.2.  3-D  graphs  for  the  optimal  solution  to  the  test  problems  with  the  second 


set  of  box  constrains 


CHAPTER  4 
CONCLUSION 


In  this  paper,  we  introduced  the  idea  of  the  active  set  strategy  from  both  the 
primal  and  dual  points  of  view  in  the  first  chapter.  Two  algorithms  are  presented  in 
Chapter  2 and  3,  namely,  the  dual  active  set  algorithm  for  linear  programming  and  the 
conjugate  gradient  active  set  algorithm  for  positive  definite  quadratic  programming 
with  box  constraints.  Some  properties  of  these  two  algorithms  are  included.  The  main 
result  is  that  both  methods  will  attain  an  optimal  solution  for  its  problem  respectively 
in  a finite  number  of  iterations.  We  implement  the  dual  active  set  algorithm  on 
some  test  problems  from  Netlib.  A comparison  of  relative  speed  of  the  conjugate 
gradient  active  set  algorithm  and  some  gradient  based  algorithms  on  a engineering 
test  problem  is  also  included. 

We  study  how  to  use  the  framework  of  the  dual  active  set  algorithm  in  Hager  [15] 
and  then  modify  it  to  solve  the  linear  programming  in  section(2.2).  Convergence  of 
this  algorithm  is  established  through  a series  of  lemmas  and  theorem  in  section  (2.3) 

The  dual  active  set  algorithm  for  the  linear  programming  problem  shares  the 
following  feature  with  the  primal-dual  method  (developed  by  Dantzig,  Ford,  and 
Fulkerson  [5]).  It  maintains  the  complementary  slackness  conditions  and  the  dual 
feasibility  while  relaxing  the  primal  feasibility  until  the  termination  of  the  algorithm. 
The  dual  active  set  algorithm  for  linear  programming  solves  the  dual  problem  of  (LP) 
starting  with  dual  feasibility  to  identify  an  active  index  set  according  to  satisfaction 
of  the  complementary  slackness  conditions.  The  dual  objective  function  value  is 
constantly  increased  by  maximizing  an  auxiliary  problem  L%  with  a set  of  active 
index  set  toward  either  the  finding  of  optimal  solution  for  the  dual  problem  and  the 
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primal  problem  (LP)  or  a conclusion  of  infeasibility  of  the  given  linear  programming 
problem  (LP). 

While  the  Simplex  method  needs  a phase  I procedure  to  find  an  initial  vertex 
to  get  started  and  the  interior  point  methods  need  a conversion  of  a general  linear 
program  into  a special  form  for  initialization,  the  advantage  of  the  dual  active  set 
algorithm  for  linear  programming  to  these  two  methods  is  that  no  further  assumption 
for  the  linear  program  is  needed.  Moreover,  any  arbitrary  initial  dual  variable  can 
be  used  to  initialize  this  algorithm. 

Although  the  dual  active  set  algorithm  possesses  similar  convergence  property  to 
that  of  the  Simplex  method  (by  the  finitely  many  possibles  of  the  different  active 
index  sets  and  the  constantly  ascent  of  the  dual  objective  function  value  with  respect 
to  different  active  index  set),  the  numerical  experiment  results  in  table  2.1  are  quite 
impressive.  That  is,  most  test  problems  have  a common  bound  on  the  ratio  of  number 
of  all  sub-iterations  and  number  of  row  (column)  of  the  constraint  matrix  A in  (LP). 
This  result  is  similar  to  that  of  the  Simplex  method.  That  is  that  the  number  of 
iterations  required  by  the  Simplex  method  is  approximately  2 to  3 times  the  number 
of  rows  of  the  contraint  matrix. 

The  conjugate  gradient  active  set  algorithm  implements  the  primal  active  set 
algorithm  presented  in  section  1.2  by  using  Conjugate  gradient  method  for  solving 
a quadratic  subproblem  with  equality  constraints  while  maintaining  an  active  index 
set.  The  algorithm  presented  in  section  3.2  is  used  to  solve  a quadratic  programming 
problem  with  box  constraints.  However,  it  is  not  restricted  to  this  type  of  problem 
only.  We  can  generalize  the  constraint  set  to  a polyhedral  type  constraints,  namely, 
{x  : Ax  < 6}. 

From  the  numerical  experiment  results  appear  in  table  3.1  through  table  3.8,  we 
find  that  the  conjugate  gradient  active  set  algorithm  is  not  much  affected  by  the  ill 
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condition  of  the  Hessian  matrix  of  the  problem.  The  performance  can  even  be  greatly 
improved  when  a modification  on  the  procedure  adding  in  more  indices  to  the  active 
index  set  based  on  the  Armijo/ Goldstein  type  line  search  routine  is  employed. 

Being  motivated  by  the  performance  of  the  conjugate  gradient  active  set  algo- 
rithm, we  look  forward  to  extending  this  algorithm  to  a more  general  function  with 
linear  constraints  and  also  to  studying  its  convergence  properties. 
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